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ABSTRACT 


Understanding turbulence degradation of electromagnetic 
wave propagation is essential for efficient operation of laser 
weapons, target designators, and imaging systems. Random 
atmospheric refractive index inhomogeneities alter the phase 
and amplitude of electromagnetic waves. 

This thesis attempts to model atmospheric turbulence 
effects by using filtered Gaussian phase screens to represent 
the random nature of refractive index changes. The simulation 
uses two-dimensional 512 x 512 fast Fourier transform (FFT) 
techniques with extended Huygens-Fresnel principles performed 
on a desk top computer. 

Simulation verification was accomplished by comparing 
calculated and theoretical spatial coherence lengths, Po. 
Phase only screens produced coherence lengths that were 30% 
larger than theoretical values. By using random phase and 
amplitude screens, the calculated coherence lengths agreed to 
within 3% of theoretical values. Saturation of the normalized 
intensity variance, c?/I?, occurred for increasing turbulence 


using a single phase-amplitude screen. 
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I. INTRODUCTION 


Turbulence in the atmosphere has a detrimental effect on 
the propagation of electromagnetic waves at optical 
wavelengths. Minute changes in the refractive index of the 
atmosphere along the path of the wave cause perturbations in 
the wavefront in both intensity and phase. Consequently, 
random atmospheric density fluctuations degrade the 
performance of optical imaging systems. Understanding this 
phenomenon is essential for future endeavors in laser weapons 
and designators, as well as imaging systems. 

This thesis attempts to model the effects of atmospheric 
turbulence using the extended Huygens-Fresnel principle. This 
theory considered both Fraunhofer and Fresnel propagation 
although the simulation focused on a single phase screen using 
the Fraunhofer approximation. 

A Rytov approximation filtered, Gaussian distributed, 
random phase and amplitude screen simulated the effects of 
turbulence on an optical beam. The magnitude of turbulence 
of this screen was proportional to the index of refraction 
structure parameter, Cn?. 

The simulation code used two-dimensional fast Fourier 
Transform (FFT) algorithms extensively to model optical wave 


propagation in the Fraunhofer regime. 


Accuracy of the simulation was checked by comparison of 
calculated and theoretical coherence lengths, Po. which were 
obtained from the e! point of the atmospheric mutual 
coherence function  (MCF). Another check verified that 
intensity variance saturation occured for an increase in 
turbulence, as measured experimentally and predicted by 


theory. 


II. BACKGROUND 


A. STATISTICAL DESCRIPTION OF ATMOSPHERIC TURBULENCE 
1. Stationarity, Homogeneity, and Isotropy 
Most theorems involving random processes require that 
the functions satisfy certain restrictions. "Stationarity" 
implies that the mean value of the random function does not 
change with time. "Homogeneity" implies that translations in 
X, y, and z coordinates (a Galilean coordinate transformation) 
do not affect variables of the random function.  "Isotropy" 
means that the random function is independent of any 
rotational coordinate transformation.  [Ref. 1] 
The last two conditions imply that the statistical 
representation of the random function for two points f (Fı) and 


f (F2) depend only on the magnitude of their separation 


r2 |. ٠ 1]‏ - وو 

If these three restrictions only hold for a local 
volume of space, L?, the random field is said to be locally 
stationary, locally homogeneous, and locally isotropic for 
| rı - T2 | < L. 

Turbulence of the atmosphere is a process that can be 
best described by a random function. Minimal violations of 
stationarity, homogeneity, and isotropy occur if the 
neighborhood is sufficiently small, less than some outer scale 


length L. [Ref. 1] 


2. The Structure Function 
Because of their nature, random functions, f, can best 
be described stochastically. The simplest moment of f is the 
mean value, «f». 
Another common parameter used to represent a random 


function is the variance, c?, defined as 


o2 = < ( f(r) = «£(r)» )? >. (2.8 


In a physical random field, such as atmospheric 
turbulence, the mean value, <f(r)>, is difficult to ascertain 
because of trends in the mean as well as variations in spatial 
position, vertically and horizontally. Therefore, the mean 
and variance are awkward parameters for describing a real, 
random, atmospheric field. To resolve this difficulty, 


Kolmogorov introduced the structure function 


De (r) = De(re-ri) = <[£(r1) - £(r2)]2>, (2.2) 


which resembles the variance at first glance. [Ref. 2] The 
structure function is an ensemble average taken over all 
possible point pairs r: and rz. [Ref. 2] 

Through dimensional analysis, and assuming isotropy, 
Kolmogorov deduced a relation between the structure function 


and the distance between the sampling points, 


Dr (r) = Cf? r2/3 (2.3) 


which holds over a limited volume, called the inertial 
subrange. [Ref. 2] Although developed for the velocity field 
v, Equation 2.3 also holds for certain atmospheric functions 
called conservative passive additives. 

Conservative passive additives are quantities that 
play no part in the turbulence of the medium. Temperature is 
not such a quantity, but potential temperature, which corrects 
for the adiabatic change in temperature with altitude, is. 


Therefore, the relation 
Dr (r) = Cr? 3ح‎ ) 2 4) 


holds within the inertial subrange where T is the potential 
temperature. 

in examining problems of optical propagation, the 
index of refraction is the essential parameter. To express 
the index of refraction as a conservative passive additive, 


it is useful to manipulate the following relation 
n = 1 + 77.6(1 + 7.52 x 10-?/^?) (P/T) x 1078 (2.5) 


where À is the wavelength of the radiation in um, P is the 
atmospheric pressure in mb, and T is the atmospheric 
temperature in K. [Ref. 3] The differential of Equation 2.5 
is 


dn = - 79 x 10-6 2*2 (2.6) 


for red light (À = 0.6 um) after ignoring the dP term 
(isobaric turbulence). Since potential temperature is a 
passive additive, the index of refraction is also. Using 
Equation 2.6 to express the index of refraction structure 
parameter in terms of the temperature structure parameter 
gives 
Das (r) » Ca? r?/3 - (79 P/T? x 10-9)? C1? r?/8 (2.7) 
[Ref. 3). 
3. The Correlation Function and Power Spectral Density 


Expanding Equation 2.2 gives 


De (r) = «f?(ri)» - 2«£(ri)f(r2)» + <£2(r2)>. (2.8) 


But if we assume an homogeneous medium, then « f?(r) » is the 
same for any point r, assuming <f(r)> = O. The structure 
function becomes 


Dr(r) - 2 Br(0) — 2 Br lr), (2.9) 


where Tatarski defined the correlation function, Br, as 


Be(r) = < £(r1) £* (r2) > (2.10) 


in a stationary, homogeneous, and isotropic medium. [Ref. 2] 

Using stochastic Fourier-Stieltjes integrals, Tatarski 
developed a relation for the three-dimensional power spectral 
density function, $(f), as the Fourier transform of the 
correlation function. Writing this in freqency (f) form as 


opposed to radian freqency (K) form preferred by Tatarski, 


a? exp (-2nif-r) B(r) d?r, (225113 


V 
and conversely, 


Bı (r) = exp (-2nif-r) é(f) d?f. (2.12) 


[Ref. 2] 
Using the fact that the correlation function is even, and 


using relations (2.7) and (2.9), the structure function for 


the Kolmogorov power spectral density is 


a (f) = 1.303 Cp? f -11/3 , (2.13) 
for the Kolmogorov power spectral density. This equation is 
analogous to 

da (KR) = 0.033 Cp? K-11/3, (2.14) 


an expression developed by Tatarski in K space that is seen 


in many publications. [Ref. 2] 


B. ELECTROMAGNETIC PROPAGATION THROUGH TURBULENCE 
1. The Wave Equation 
Central to any study of electromagnetic propagation 


are the four equations of Maxwell, presented here in Gaussian 


units 


(2215) 


4 
اب‎ 
Il 
O 
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-ikn2E, (2.16) 


V + (n2E) = 0, (2.398 


V x E = ikH, (2.18) 


for the assumptions that the medium has zero conductivity, 
unit magnetic permeability, and sinusoidal dependence. 


Taking the curl of Equation 2.18 gives 
V x (V x E) = V x (ikH). (2.19) 


Using a vector identity and Equation 2.16, Equation 2.19 
becomes 


-V2E + V(V + E) = k2n2E. (2.20) 


Expanding Equation 2.17 gives a relation for 
V + E, which when inserted into Equation 2.20 gives the wave 
equation 


V?E * k?n?E + 2V[E + Vílog n)] = 0. (2.210 


If the wavelength of the electromagnetic wave is small 
compared to the dimensions of the refractive inhomogeneities, 
then the 2V[E - V(log n)] term is negligible. In this case, 


the wave equation reduces to the simple form 


V?E + k2n2E = 0. (2.22) 
[Ref. 3] 
2. The Born Approximation 
One method of solving the wave equation, employed by 
both Tatarski [Ref. 2] and Clifford (Ref. 3], is the method 


of small perturbations or the Born approximation. This method 


expands the electric field and index of refraction as a power 
series 
E = Eo + Ex: + °°°, (2.23) 


ho- Jot nit o, (2.24) 


where Eo is a constant in time, and E, and nm: are time 
varying. 
Inserting these two relations into Equation 2.22 and 


equating same order terms gives 


V?Eo * k?Eo = O, (2.25) 


V?2E1 + K?E1 + 2nik*Eo = O, (2.26) 


ignoring all second and higher order terms. 
Assuming that the electromagnetic wave propagates in 


the z-direction and E o = exp(ikz), then Equation 2.26 becomes 
V?Ei * k?E: -* 2nik?exp(ikz) = 0. (2 27) 


Clifford solved this wave equation for E: by utilizing 
a Green's function. The solution is the convolution of a 
plane wave Green's function with the source term, shown here 
in scalar form 
1 exp[ik|r-r' |] 


0 A A o AO E ep Ik :8ه(‎ 
41 V FE 


Eı (F) 


(2528) 


where r is the position of a specific point in the image plane 


and r' is the postion of specific source point. [Ref. 3] 


Figure (2.1) shows the coordinate system used to express the 
solution to the wave equation for propagation between planes 
Z' and P' to z and P- 


For normal propagation situations, |z-z'| = R >> | م‎ - 


so that Equation 2.28 expands into the form‏ ]|'م 


1 exp{ikR (1 + | ې ۶۰ 'م -م‎ ... [( 
Eilr) 3 ۳ — = 
4n Jv R (1 + |p-p'|?/2R? + ce.) 
X [2k2 ni(r') exp(ikz')] ar, (2.29) 


using the binomial expansion. Dropping those terms consistent 
with the Fresnel approximation in Huygens-Fresnel optics 


theory, this equation reduces to 


k2 exp (ikz) 
EME). = _ لر‎ exp{ik| p-p ' |?/2R} ni(r') d?r, 


AnR ۷ 
(2.30) 
which is known as the Fresnel diffraction formula. (Ref. 3] 
3. The Rytov Approximation 
Another method of solving the wave equation, also used 
by botk Tatarski and Clifford, is the method of smooth 
perturbations or the Rytov approximation. This technique 


assumes that the electric field is of the form 


E = exp(¥) = exp(X+iS) = A exp(iS). )2.31( 
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Image Plane 


Aperture Plane 





Figure 2.1 Aperture and Image Plane Geometry. 
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Using power series expansions for the electric field, 
Equations 2.25 and 2.26 still hold. Placing Equation 2.31 
into these two gives 


(V2E)/E * k*n*(r) = 0. (2.32) 


Tatarski shows that this method gives the same results as the 
method of small perturbations but the range of validity is 
larger than for the Born approximation. Equation 2.26 holds 


as long as V-E is small compared to A. (Ref. 2] 


C. HUYGENS-FRESNEL THEORY 

Another way of obtaining a relation for the electric field 
amplitude at the observation plane, is by extending the ideas 
of the Huygens-Fresnel theory to include a random phase tern. 
The Huygens-Fresnel theory offers this solution for the 


electric field seen in the observation plane 


E(r) = - ik E(r') exp (ik|r-r' |) do '. (2.399 


2n S Ir-r'| 
(Ref. 5] 

Lutomirski and Yura develop an extension of the Huygens- 
Fresnel idea which incorporates a random phase term which is 
equivalent to the form used in the Rytov approximation, 
exp(¥), where Y isa complex phase. [Ref. 4] This extended 


integral is 


E(F) = - ik E(r') exp( V (r)) exp(ik|r-r'|) dp '. 
2n 5 |Ir-r'| 
(2.34) 
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In the geometrical optics limit, V(r) becomes the optical 


phase 


V (r) = ik | m(Zz) dz. (2.35) 


Expanding Y into a power series, Equation 2.35 reduces to 


Equation 2.30 as developed by Tatarski and Clifford. 


D. MUTUAL COHERENCE FUNCTION 

In trying to get a quantitative measure for turbulence, 
Lutomirski and Yura [Ref. 4] examined the average intensity 
of the wave at the observation point. The average intensity 
at r is defined as 


<I(r)> - « E(ri1) E*(r2) >. (2.36) 


From Equation 2.34: 


<I(r)> = < (k2/4n2) E(ri')E*(rz') exp( Y (X1) + Y* (F2)] 


X expl ik(]|r:-r; ' | = |rz-rz'|)] dos dp 2 > 
[ri-ra] Era 
)2.37( 
Picking out only the time-dependent portion of this result, 


the entire integral reduces to evaluating 
« exp( V (r1) + V*(r:)] ». (2.38) 


This term is the atmospheric mutual coherence function (MCF) 
which contains the elements of the propagation through 
turbulence. 

Lutomirski and Yura (Ref. 4] relate the atmospheric MCF 


to the wave structure function, D, by the relation 
13 


MCF(p) = «exp[ V (rji) + W*(F2)]> = exp[-D(o )/2) (2.39) 


where D(o) = Dx (0) + Ds (p) from the Rytov relation Y = X+iS 


or 
MCF ( 0 ) = expl-( 0 / po)?/?], (2.40) 
where 
po = [ 1.46 k2 R Cn? dz )-3/5, (2.41) 
0 
and k is the wavenumber. Cn? is the refractive index 


structure parameter defined by Equation 2.7, and R is the 
distance of propagation through the atmosphere. The coherence 
length, Po. is defined as the distance where transverse 


spatial coherence of the wave drops to the e^i. [Ref. 1] 


E. THE DIFFRACTION INTEGRAL 
Like Lutomirski and Yura, Roberts [Ref. 6] begins with the 


Huygens-Fresnel integral for the propagation of light waves 


after making the Fresnel approximation انا‎ i = R >> iz n 
E(r) - -ik do E(r') exp{ik[R?+| p =p [2] 972 eee 
AnR S 


In this expression, k is the wavenumber and the notation is 
the same as shown in Figure (2.1). 
Employing the standard practice of expanding by the 


binomial expansion and neglecting smaller terms gives 


E(Y) = -ik exp[-ikR] dp ' E(r') عدا صعدء‎ Pp E 
2017 5 2R 
(2.43) 
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Since the exponential term outside the integral doesn't affect 
intensity measurements, it can be neglected. Expanding the 


quadratic term gives the relation 


E (r) EE dp ' E(r') exp{ik 0٦ on 
2R 


2018 5 





-ik explik 0'2] do' E(r') explik 2] expl-ik0O - D'J. 
21718 m a 2R R P 


(2.44) 
This is called the convolution form of the Fresnel diffraction 
integral which differs from the Fraunhofer approximation only 


by the quadratic phase term exp[ikpo?]. [Ref. 6] 
2 
If we denote the Fourier transform of E(r) as E(f), then 


E(f) = | dp exp(-2nif-6 ) Elx). (Pas) 


5 


Inserting E(r) as defined in Equation 2.43, gives 


E(f) = -ik |dp' E(r') |dp exp(-2nif-6 ) explik| 6 '- 6 |?- 
2nR Js 5 AR 

(2.46) 

Changing variables to q = p - م‎ , the electric field 


spectrum becomes 


E(f) = -ik dp' E(r') exp(-2nif:f) 
2nR J, f P 
x [as exp (-2nif-q) exp(ik q?). (2.495) 
S AR 


The q integral is shown to be the Fourier transform of the 


Gaussian function 


و 


2niR exp(-2n?iR£?). (2.48) 
k k 





(Ref. 6] Substituting this function into Equation 2.47 gives 


E(f) = exp(-2n? iR £?) do ' E(r') exp(-2nif-P), (2.49) 
k 


which is the inverse transform of what we sought, 
E(r) = fas exp(2nif-r) exp(-2n2iR f ) 
k 


X fa " EU exp (-2nif.p ) (2.50) 


This is the transfer form of the Fresnel diffraction integral. 
[Ref. 6] 

Equations 2.44 and 2.50 are two different forms of the 
diffraction integral that are useful for two types of wave 
propagation. The convolution form, Equation 2.44, is best 
suited for long distance propagation (Fraunhofer regime) since 
R is in the denominator of the exponential term. The transfer 
form, Equation 2.50, is best suited for short distance 
propagation since R is present in the numerator of the 


exponential term. [Ref. 6] 
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III. NUMERICAL SIMULATION 


A. SIMULATION HARDWARE AND SOFTWARE SPECIFICS 

This numerical simulation modeled the propagation of an 
electromagnetic wave through a random turbulent medium. The 
procedures for this simulation required the creation and 
manipulation of multiple two-dimensional, 512 by 512, single 
precision, phase screens, requiring a minimum memory of 
several megabytes. 

A COMPAQ Deskpro 80386-20 computer configured with a 16 
megabyte RAM and 64 megabyte hard drive met hardware 
requirements. The computer also had EGA-VGA graphics and both 
HP Laser Jet II and Panasonic KX-P1092i printers. A Weitek 
1167 math coprocessor increased execution speeds of the 20 Mhz 
80387 coprocessor by a factor of 3-4. 

The software requirements for writing the simulation code 
were met with the MicroWay 32 bit NDP FORTRAN-386 compiler 
and the Phar Lap operating system extender. It was chosen 
over the Science Applications International Corporation's SVS 
Fortran-386 compiler which had 10-20% faster program execution 
times but numerous bugs present in the graphics routines. 

Unfortunately, due to limitations of version 1.4e of the 
NDP Fortran-386 compiler, the full graphic capabilities of VGA 


were not supported, so EGA graphics were used for screen 
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displays. Hardcopy graphs and plots were accomplished using 


Grapher and Surfer programs from Golden Software. 


B. REQUIRED ALGORITHMS 

For the successful operation of such a model both the 
random number generator and the two-dimensional fast Fourier 
transform (FFT) routines needed speed and accuracy. These 
will be discussed in detail in the following sections. 

1. Random number generator 

The effects of turbulence on the propagation of an 
electromagnetic wave are of a random nature. To emulate this 
random nature required a pseudo-random number generator 
algorithm. The minimum criteria for a quality random number 
generator was that the products be: (1) sufficiently random, 
(2) reproducible, and (3) rapidly obtained. 

The NDP FORTRAN compiler did not come equipped with a 
random number generator, so an algorithm of the linear 
congruent type especially designed for 32 bit machines was 
implemented as a subroutine attached to the main program. 
This had the benefit of not binding this portion of the code 
to a specific machine setup. This random number generator 
algorithm came from Dr. Milne, who obtained it from class 
notes prepared by Dr. Harrison [Ref. 7]. This algorithm 
appears as subroutine RAN (Appendix A) in the simulation code. 

The first criteria, that the generator produce numbers 


sufficiently random, was hard to quantify. A series of tests 
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exercised the statistical independence of succeeding numbers. 
There are many such tests that exist, the most common of which 
is the "Chi-Square" test. The following is a list of all 
tests used to check the quality of the random number generator 
algorithm: 


1. Frequency Test. This tests the uniformity of the 
sequence of numbers. 


2. Serial Test. This tests the two-dimensional uniformity 
of the sequence of numbers. 


3. Runs Up and Down Test. This tests to see how long 
sequences of random numbers are that either go 
successively up or down. 

4. Lagged-Product Test. This tests for correlations 
between random numbers Ri and Ri:3 where Jj is the given 
lag value. A lag of 3 is especially difficult to pass. 

5. Repeat Test. This simple qualitative check is to 
determine how long the sequence of random numbers 
becomes before the first number is repeated. 

6. Scatter Plot Test. This is another qualitative test 
where random numbers plotted on a x-y plane are visually 
checked for correlation. 

These six tests were performed for a variety of seed 
values to determine which gave the best performance. Results 
of the first five tests for two good seeds can be found in 
Table 3.1. 

Choosing the same initial seed value gave 


reproducibility of the random numbers. The algorithm's simple 


two lines of code assured rapid execution. 
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TABLE 3.1 RANDOM NUMBER GENERATOR TEST DATA 


This table provides the results of statistical tests of 
the random number generator algorithm. The X? values came 
from the Chi-Square table in Bevington (Ref. 11] 

The results of the Lagued Product test correspond to the 
theoretical mean, ur, the calculated mean, u, the theoretical 
standard deviation, or, and the calculated standard deviation. 


Degrees of 
Test Freedom Probability 
“45739. 91 % 
Frequency 
94377 اتد‎ 82 # 


45739 








Serial 
94377 


- 45739 | 0.250 | 0.230 0.088 0.092 
2. Fast Fourier Transform 
Since the most often used algorithm in the simulation 
was the Fast Fourier Transform (FFT), efficiency was critical. 
Two different routines were compared. 
The first FFT routine considered was a package offered 
by NDP of one- and two-dimensional machine language FFT 
routines designed for easy linking by the NDP FORTRAN 


compiler. This package had routines available for real and 


complex arrays. 
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The second FFT routine considered was a simple real 
array subroutine coded by Dr. Walters from a BASIC 
demonstration package provided by Infotek. 

Each FFT routine was compared for numerical output and 
speed of operation for the one-dimensional case. As expected, 
all routines were virtually identical in numerical output. 
However, speed of operation differed widely, as seen in Table 


2. 


TABLE 3.2 FFT SPEED OF OPERATION COMPARISON 


This table contains results of FFT execution time 
comparison between NDP machine language routines and Dr. 
Walters' subroutine. Since NDP FORTRAN does not have a timing 
function with sufficient accuracy, execution times listed are 
mean times found from several runs, timed by stopwatch. 


NDP Routines Walters' Routine 
CFFT RFFT w/o € w/ Weitek 


Em ———‏ ہے۔ : aaa‏ لي ——— وس س سم A‏ 


e (s) 
2.84 





The differences in run time for the two NDP machine 
language routines were due to the CFFT routine being more 
cumbersome in converting repeatingly between real and complex 
arrays. The compiled subroutine was significantly faster than 
both CFFT and RFFT NDP machine language routines unless the 


Weitek was deactivated. Using a subroutine enclosed within 
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the main program had the additional benefit of not being 
machine dependent. Therefore, it was decided that Dr. 
Walters' subroutine would be used for all FFTs needed in the 
simulation code. Conversion of this subroutine for two- 
dimensional use was relatively easy. The NDP two-dimensional 
FFT routines never worked correctly. 

A FFT "breaks down" a function into its cosine and 
sine constituents (real and imaginary terms, respectively). 
The range of spatial frequencies represented in a FFT depend 
on the array size chosen. The minimum frequency is fmin = 1/W 
and the maximum frequency is fmax = n/W, where W is the 
physical width of the array and n is the number of sampling 
points across the array. When energy at frequencies near 
these cutoff values appeared, problems developed due to the 
discrete nature of the FFT. 

One problem exists in representing spatial frequencies 
smaller than fmin due to the high amplitude, low frequency 
"tilt" of the f -!!/* filtered wavefront. Since a tilted 
straight line approximates a cosine or sine function for only 
a small region, it helps to make the wavefront smaller than 
the FFT array size. Hence, to help represent lower 
frequencies in the simulation coding accurately, an aperture 
array was a user chosen variable smaller than the FFT array 
size. 

Another problem that manifests itself is "aliasing", 


which occurs when high frequency features of a function are 
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missed by having too few sample points. This means that 
important information for the function present in the 
frequencies greater than feax are folded back and appear at 
lower frequencies. Choosing the distance between sample 
points smaller than the coherence length of the propagating 
electromagnetic wave, found theoretically by Equation 2.39, 


alleviated this problen. 


C. SIMULATION DESCRIPTION 
1. User defined quantities 
For operation of the simulation program, the user had 
the option of varying several quantities directly from the 
keyboard: 


1. Array size. This can be varied up to a maximum of 512 
by 512 points. 


2. Aperture sampling points. This can be varied up to a 
maximum of 512 by 512 points. Ideally it should be 


smaller than the array size to avoid "tilt" effects, but 
large enough to avoid “aliasing” effects. A size one- 
fourth to one-half the array size works best. 


3. Aperture shape. A choice of square or circular aperture 
was provided. 


4. Seed. The pseudo-random number generator required a 
five digit input seed. Table 3.1 lists two quality 
seeds. All data was obtained using the seed 94377. 


5. Turbulence parameter. Any value of Cn? above zero can 
be chosen. Typical values of Cn? were 10-17 through 
10-13 m-2/3. 


6. Wavelength. Any value of wavelength for the electro- 
magnetic wave can be chosen. A value of 0.6 yum was used 


throughout. 
7. Distance. This was the distance from the aperture to 
the observer. This simulation used 1 km. 
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8. Aperture width. Adjustable to ensure that "aliasing" 
effects do not occur for a specific choice of Cn?. 
Typically, a 100 mm aperture was chosen with a path 
length of 1 km. 

For the accumulation of multiple data points for 
statistical purposes, a slightly different version of the 
program read user choices from a disk file. 

2. Aperture Mutual Coherence Function 
Computing far-field diffraction patterns of intensity 
for an electromagnetic wave through the aperture with no 
turbulence gave the Mutual Coherence Function (MCF) of the 
aperture. Since we assumed a point source at infinity, the 
electric field wave front was planar at the aperture. The 
electric field that passes through had a value of one within 
the boundaries of the aperture, and a value of zero outside 
the boundaries. The FFT subroutine computed the diffraction 
pattern, the electric field present at the image plane. The 
conjugate square of this electric field was the intensity 
field that would be seen by an observer. Taking the inverse 
FFT of the intensity field gave the aperture MCF, which was 
identical to the autocorrelation of the aperture function. 
3. Phase screen generation 
Simulation of effects of turbulence requires that the 
aperture planar electric field be multiplied by the phase 


screen. Two methods exist to create the phase screens. 
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a. Classical method 

In this method n? Gaussian distributed random 
numbers filled a n by n real array. Since the pseudo-random 
number generator RAN produced uniformly distributed numbers, 
a subroutine GAUSS transformed them to a Gaussian distribution 
with unit variance, Knuth [Ref. 8]. This subroutine appears 
in Appendix B. Taking a direct FFT creates two Gaussian 
distributed arrays, phaser and phasei, which are the real and 
imaginary spectral amplitudes of the original random array. 
Since these two arrays are in spatial frequency space, 
multiplication by Equation 2.13 occurs. An inverse FFT 
returns the filter to real space. Because of the symmetry of 
the filter, discussed later, only one array in real space 
results, containing all n? pieces of information produced by 
the subroutine RAN. 

Martin and Flatté method‏ ط 

In this second method suggested by Martin and 
Flatté (Ref. 9], two n by n arrays, called phaser and phasei, 
are each filled with n? Gaussian distributed random numbers. 
These are considered the real and imaginary components in 
Spatial frequency space. Filtering occurs as f-!!/3 then an 
inverse FFT returns the array to real space. Because no 
information is lost in filtering, two arrays result, each with 
n? pieces of information. Only one of these arrays need be 
used for the following simulation code. This method requires 


twice the random numbers as the previous method but GAUSS 
\ 
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produces two random numbers and half as many two-dimensional 
FFT steps are required. For use in Fresnel propagation, where 
multiple phase screens exist, this method is certainly the 
most efficient. 
4. Filtering 

Phase screens need to be filtered in frequency space 
to ensure the proper power law form. From Tatarski [Ref. 1] 
and Martin and Flatté [Ref. 9], the correct form for the 


filtering function, 4:, is 


Pr = 2n k? öz On, )3 


where k is the wavenumber of the electromagnetic wave, 6$; is 
the slab thickness, and @ is the power spectral density 
function. Using Equation 2.13 for the power spectral density 


function and the definition of wavenumber, k = 2n/À, gives 


Pi = 1.303 223 ٨٨۸/٢٤٢ CHE o 7 )3-2 


which is the filtering function used in the subroutine FILTER 
in the simulation code. Arrays phaser and phasel are both 
multiplied by the square root of Equation 3.2 to obtain a 
factor for the real and imaginary amplitudes so that the 
intensity has the proper f-!!/? power law. 

Although this filtering scheme looks simple, a few 
subtle points make the process much more complicated than is 
immediately evident. In fact, perfecting this filtering 


process constituted the major effort of this thesis. 
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First of all, most of the relations developed in the 
background references were in terms of the radian frequency 
K. The FFT routine used was in terms of the spatial frequency 
f - K/2n. Consequently, all relevant equations had to be 
expressed in this latter form, requiring the tracking down of 
many 2n terms. 

Secondly, due to the two-dimensional isotropy of 
turbulence, the correct filtering scheme is of a circular 
nature. This means that the frequency seen in Equation 3.2 
is really f = (fx2 + f,2)1/2, which is radial everywhere but at 
the zero frequency (DC term) where it is has a value of zero 
[Ref. 9]. 

For convenience of coding, the zero frequency should 
be at the center of the arrays phaser and phasei. But in the 
FFT algorithm, spatial frequencies are arranged in a manner 
that is quite different. The zero frequency term is at the 
upper left corner of the array. The Nyquist "folding" 
frequency is present as a cross that divides the arrays into 
four sections. So, prior to filtering, a shuffling of data 
points is required as discussed by Brigham [Ref. 10] and as 
seen in Figure 3.1. After filtering, an inverse shuffle 
returns all frequencies to their previous locations. 

Another important consideration in the filtering 
process was that frequency units be correct according to the 
physical dimensions of the aperture. This required that the 


frequency array have a normalization factor of 1/NA where N 
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is the total number of points across the array, and A is the 
sampling interval (physical units between the aperture array 
points). 
5. Propagation 

The Huygens-Fresnel techniques formed the basis for 
the simulation's treatment of optical propagation. For 
simplicity this simulation only considered Fraunhofer 
propagation. The addition of subroutines to handle the 
dynamics of Fresnel propagation should not be too difficult 
other than the need for dynamic scaling of the numerical mesh 
to account for diffraction spreading. 

Martin and Flatté provided a procedure to simulate the 
propagation of a wave in a turbulent medium. [Ref. 9] 

1. Add the random turbulence effects to the electric field 


distribution by meshing the aperture electric field 
function with a random phase term exp(?): 


E(T') > E(£') exp(V ). (3.3) 

2. Take the Fourier transform of the result of Equation 3.3 
to obtain 

E(f') = DFTIE(F')). (3.4) 


3. Add a term to the result of Equation 3.4 to obtain 


E(F) = E(f') exp{-2niRf |} (3.5) 
k 


which is the exact same relationship as in Equation 
2.48 that was described by Roberts as the 

convolution form of the Fresnel diffraction integral. 
[Ref. 6] 


4. Take the inverse Fourier transform of Equation 


3.5 to get an expression for the electric field 
distribution at the image plane: 
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E(r) = IFT{E(f)} (3.6) 


If one were taking a multiple screen approach to 
solving for the propagation of the wave as suggested by Martin 
and Flatté, the entire process would be reiterated for as many 
segments of length R as required. In this numerical 
Simulation, only one segment was used for all propagations. 
Adding this iterative section to the code would be a 
straightforward modification if dynamic array sizes were 
incorporated. [Ref. 6] 

For distances of propagation where the phase front 
cannot be assumed planar, an additional quadratic term would 
have to be added at step 3 for the transfer form of the 
Fresnel diffraction integral, a= Seem in Equation 275 

6. Atmospheric Mutual Coherence Function 

We computed the atmospheric mutual coherence function 
after meshing the random phase screen with the aperture 
electric field function (after step 1l ٢ ٢٢٢ 5 and Flatté 
method). As with the aperture MCF, the next step was 
computing the diffraction pattern using the FFT routine, which 
represents the complex electric field that an observer would 
see because of the aperture and the turbulence. The effects 
of turbulence tends to spread out the diffraction pattern so 
that resolution would be lost for any optical imaging system. 
Taking the conjugate square of this diffraction pattern gives 


the intensity field, which negates the additional exponential 
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term present in Equation 3.5. Taking the inverse FFT of the 
intensity field gave the total mutual coherence function. 
Dividing the total MCF by the aperture MCF leaves the 
atmospheric MCF. 
7. Coherence Length 

As stated previously, the coherence length, Do, was 
the distance transverse to the direction of propagation where 
the spatial coherence of the «E*E» wave dropped to e-!. 

For a continuous plot of the atmospheric MCF, the 
e”! distance could be picked off the curve easily. Due to 
the discrete nature of the simulation, the MCF curve was not 
continuous, but a series of points equal to the choice of the 
pixel width of the aperture array. Therefore this simulation 
used a linear interpolation routine to pick out values for 
the coherence length as seen in the results of the following 


section. 


IV. RESULTS 


In order to verify that this numerical simulation produced 
correct results, we had to find ways to check the validity of 
the code. in the code, presented in Appendix C, there are 
several checkpoints where the accuracy of the simulation could 
be verified. The following sections discuss these checkpoints 


in increasing order of sophistication. 


A. APERTURE MCF 

An initial point to check the accuracy of the simulation 
code was to look at the mutual coherence function (MCF) of 
the aperture. Since the code offered two choices of aperture 
shape, square and circular, both were verified. 

Using the two-dimensional FFT routine, the diffraction 
patterns for 100 mm wide square and circular apertures were 
calculated. The square aperture diffraction pattern in Figure 
4.1a displays the familiar [(sin x)/x)? form discussed by 
Hecht with the correct maximum value and node spacing. (Ref. 
5] The circular aperture diffraction pattern in Figure 4.1b 
displays the Airy pattern derived from the J: Bessel function 
solution carried out by Hecht.  [Ref. 5) 


An inverse 2-D FFT of each 100 mm wide aperture shape 


produced the aperture MCF (identical to aperture 
autocorrelation). The square aperture autocorrelation in 
Figure 4.2a 15 the predicted triangle function. The 
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Figure 4.2 Autocorrelation (Aperture Mutual Coherence 
Function) Of A 100 mm Wide (a) Square Aperture And A 
100 mm Wide (b) Circular Aperture 
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circular aperture autocorrelation in Figure 4.2b is as 
expected for a "top hat" function. 
The algorithm gave the correct MCF results which verified 


the simulation up through the aperture MCF point. 


B. FILTERING 

The most crucial portion of the entire simulation code was 
the correct frequency filtering algorithm. A succinct way of 
verifying this came from examining the result of taking the 
logarithm of both sides of Equation 3.3 

log{®r} = log{1.303(2n)3%(1/\)2Cn28:f -11/3} 
which is 

logit} = 1log{1.303(2n)$(1/X\)?Cn28:} - (11/3)loglfl, (4.1) 
an equation for a straight line with slope of -11/3. 

Figure 4.3 shows a plot of filtering function vs. 
frequency f with a value 10-!? for Cn? and a wavelength of 600 
nm (red light) using the Martin and Flatté method of filtering 
discussed earlier. A linear fit to the data points down to 
the Nyquist frequency had a slope of -3.88. This was a 
deviation from the predicted value of -11/3 (-3.83) by 1.3%. 
Therefore, it appears that the filtering is being accomplished 
correctly. The upturn of points occurs at the Nyquist folding 


point that is present in the array scheme for the FFT routine. 


C. ATMOSPHERIC MCF 
An important place to check results of the simulation was 


at the calculation of the atmospheric MCF. Following the 
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Figure 4.3 Slice Through Phase Screen Function 
Filtered By The Martin And Flatté Method Showing 
Its Dependence On Spatial Frequency. 
Slope Of -11/3 Comes Directly From The Kolmogorov 
Power Law. The Upturn At Higher Frequencies Is Due 
To Nyquist Folding Present In The Phase Screen Array Code. 
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procedure for calculating the atmospheric MCF outlined in 
section III, first the diffraction pattern was calculated. 
Figure 4.4a shows the diffraction pattern for a 100 mm wide 
square aperture that results for red light (À = 0.6 um) with 
Cn? = 10-16. Notice that the diffraction is only slightly 
affected from what is shown in Figure 4.1a. As the turbulence 
structure parameter Cn? increased, the diffraction pattern 
began to spread out more until it no longer resembles the 
[(sin x)/x]? form it had without turbulence, as seen in 
Figures 4.4b, 4.4c, and 4.4d for values of Cn? = 10-15, Cn? = 
10-14, and Cn? = 10-!3 respectively. Figure 4.4c shows that 
modest amounts of turbulence displace the image centroid, due 
to low frequency tilting of the electric field. Since the 
diffraction pattern is analogous to the intensity that would 
be seen at the image plane, this clearly shows the blurring 
Of an image commonly seen looking through a turbulent 
atmosphere as well as the "dancing" image effect at lower 
levels of turbulence. 

The atmospheric MCF was calculated for different values 
of turbulence by the method described in section III using the 
classical method of filtering. These are seen in Figures 4.5 
for Cn? values of 0, 10-16, 10-15, 10-14, and 10-182. With no 
turbulence, the atmospheric MCF is a value of 1.0 everywhere, 
seen as a straight line. With increasing turbulence, the 


atmospheric MCF curves fall off rapidly to zero. 
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(e) (d) 


Figure 4.4 Diffraction Patterns For A 100 mm Wide Square 
Aperture For Four Different Values Of Turbulence: 
(a) Ca?210-!$ (b) Cas?210-!5 (c) Ca,?-210-!* (d) Ca?z10-!3., 
Displacement Of The Beam Centroid Is Obvious In (c) Due 
To Low Frequency Tilting Of The Electric Field. 
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Figure 4.5 Atmospheric Mutual Coherence Function Curves 
For Five Different Levels Of Turbulence. 
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D. COHERENCE LENGTH 

The e-! point of atmospheric MCF curves give the coherence 
length, po. In the simulation, we achieved this by a method 
of linear interpolation as discussed previously. One should 
not put too much faith in a statistical sample of one, 
however. For this reason, a modification of the simulation 
was made to run several samples of each different set of input 
parameters so that a better statistical representation of 
results occurred. 

Theoretical coherence lengths were calculated from 
Equation 2.41 [Ref. 1] for three different values of Cn? using 
R-1000 m and A-0.6 um. Several runs of data were produced to 
investigate the effect of different sizes of both the electric 
field array and the aperture array. Calculated results are 
discussed and displayed in following pages. 

Figure 4.6 shows the trend for coherence length values as 
a function of aperture array sample points using the classical 
method of filtering for Cn? = 10-!3. Ten samples were obtained 
for each data point. As discussed in section III, aliasing 
makes the values of po too high if there are too few sample 
points. Using aperture arrays no larger than one-half the 
width of the main array was done to avoid tilt problems. The 
curves asymptotically approach coherence length values of 4.18 
mm for the 512 x 512 array and 4.49 mm for the 256 x 256 
array. These represent errors of 39.3% and 49.7% 


respectively. 
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Figure 4.6 Calculated Coherence Length Values For 
Cn? = 10-13 Using Different Aperture Array Sizes 
And Classical Method Of Filtering. 

Error Bars Shown Are For Case Of 512 By 512 Array Size, 
Representing The Standard Deviation For Ten Samples. 
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Figure 4.7 shows coherence length trends for Cn? = 
10-14 also using the classical method of filtering. Ten 
samples were obtained for each data point. Results were 15.78 
mm for the 512 x 512 array and 18.43 mm for the 256 x 256 
array, representing errors of 32.1% and 54.2% respectively. 

Figure 4.8 shows coherence length trends for Cn? = 
10-15 also using the classical method of filtering. Ten 
samples were obtained for each data point. Results were 61.99 
mm for the 512 x 512 array and 76.27 mm for the 256 x 256 
array, representing errors of 23.3% and 60.3% respectively. 

Figure 4.9 shows coherence length trends for Cn? = 10-13 
using the Martin and Flatté method of filtering. Ten samples 
were obtained for each data point. Results were 4.07 mm for 
the 512 x 512 array and 4.46 mm for the 256 x 256 array, 
representing errors of 35.7% and 48.7% respectively. 

Figure 4.10 shows coherence length trends for Cn? = 10714 
also using the Martin and Flatté method of filtering. Ten 
samples were obtained for each data point. Results were 14.94 
mm for the 512 x 512 array and 20.46 mm for the 256 x 256 
array, representing errors of 25.0% and 71.2% respectively. 

Figure 4.11 shows coherence length trends for Cn? = 10715 
also using the Martin and Flatté method of filtering. Ten 
samples were obtained for each data point. Results were 60.74 
mm for the 512 x 512 array and 74.04 mm for the 256 x 256 


array, representing errors of 27.7% and 55.6% respectively. 
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gure 4.7 Calculated Coherence Length Values For 
Cn? = 10-14 Using Different Aperture Array Sizes 
And Classical Method Of Filtering. 
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Error Bars Shown Are For Case Of 512 By 512 Array Size, 
Representing The Standard Deviation For Ten Samples. 
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Figure 4.8 Calculated Coherence Length Values For 
Cn? = 10-15 Using Different Aperture Array S1zes 
And Classical Method Of Filtering. | 
Error Bars Shown Are For Case Of 512 By 512 Array Size, 
Representing The Standard Deviation For Ten Samples. 
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Figure 4.9 Calculated Coherence Length Values For 
Cn? = 10-1% Using Different Aperture Array Sizes 
And Martin € Flatté Method Of Filtering. 

Error Bars Shown Are For Case Of 512 By 512 Array Size, 
Representing The Standard Deviation For Ten Samples. 
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Figure 4.10 Calculated Coherence Length Values For 
Cn? = 10-14 Using Different Aperture Array Sizes 
And Martin & Flatté Method Of Filtering. 
Error Bars Shown Are For Case Of 512 By 512 Array Size, 
Representing The Standard Deviation For Ten Samples. 
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Figure 4.11 Calculated Coherence Length Values For 
Cn? = 10-15 Using Different Aperture Array Sizes 
And Martin € Flatté Method Of Filtering. 

Error Bars Shown Are For Case Of 512 By 512 Array Size, 
Representing The Standard Deviation For Ten Samples. 
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The most accurate results were obtained when the aperture 
array Was one-fourth the size of the electric field array, as 
seen in Figures 4.6 through 4.11. For this reason all 
subsequent coherence lengths were calculated using a 512 by 
512 array with a 128 by 128 aperture array. 

Figure 4.12 shows coherence length values calculated with 
different values of turbulence using the classical method of 
filtering. Ten samples were obtained for each calculated 
value. Although the calculated coherence length values 
deviated from theoretical values by approximately 30%, a 
linear fit to the data shows that the slopes were correct, 
implying a constant error in the algorithm due to 
underestimating turbulence effects. 

Adding a log-normal, random-amplitude screen to the GAUSS 
subroutine was attempted to improve the simulation results for 
coherence lengths. Having both phase and amplitude screens 
present is equivalent to the complete Rytov approximation, 
Equation 2.31. Eight samples were obtained for each 
calculated value. Figure 4.13 dramatically shows that this 
simple addition produced calculated coherence lengths that 
were within 3% of theoretical values, from a linear fit to 
the data. 

Table 4.1 presents coherence lengths calculated by the 
simulation using both random phase and amplitude screens in 


the complete Rytov approximation. All data points were 


48 


100 


سم 
© 





Coherence Length (mm) 


es Calculated coherence length » 
Theoretical coherence length D 





Figure 4.12 Calculated Coherence Length Values As A 
Function Of Cn? With Only Random Phase 
Screen In The Simulation. 
Calculated Values Are About 30% High For All Levels Of 
Turbulence, From Linear Fit To The Data. 
Error Bars Represent Standard Deviation For Ten Samples. 
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Figure 4.13 Calculated Coherence Length Values As A 
Function of Cn? With Both Random Phase And Amplitude 
Screens In The Simulation. 
Calculated Values Are About 3% High For All Levels Of 
Turbulence, From Linear Fit To The Data. 
Error Bars Represent Standard Deviation For Eight Samples. 
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obtained using a 512 by 512 array and 128 by 128 aperture 


array. 


TABLE 4.1  COHERENCE LENGTHS 


This table shows theoretical values of coherence lengths 
calculated from Equation 2.39 for eight different values of 


Cn?. 
theoretical Po (mm) calculated po (mm) 
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E. INTENSITY VARIANCE SATURATION 

The last checkpoint available for verifying proper 
operation of the simulation was whether the code led to the 
saturation of intensity variance for increasing turbulence. 
This should happen because the Rytov assumption of an 
exponential random term, exp(¥), has a magnitude bounded by 
plus or minus one. Figure 4.14 shows that saturation does 
occur around Cn? = 10-!4 at a value of unity, as expected. The 
normalized variance > 1 present at this level of turbulence 
is due to a "dancing" beam centroid caused by low frequency 
tilt of the electric field as discussed previously. A decline 
to unity normalized variance appears in the calculations for 


higher turbulence values. (Ref. 12] 
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Figure 4.14 Normalized Intensity Variance As 
A Function Of Cn?. 
Saturation Occurs Around Cn? = 10-14. 
Large Variances Near Saturation Are Due To Beam 
Centroid Displacement Cause By Low Frequency Tilt. 
Error Bars Represent Standard Deviation For Ten Samples. 
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V. CONCLUSIONS 


The program coding discussed in this thesis simulated the 
propagation of optical waves through a random media using two- 
dimensional fast Fourier transform (FFT) techniques. The 
coding used Fraunhofer propagation throughout. An extension 
to include multi-step Fresnel propagation is straightforward. 

Aliasing problems in the FFT were avoided by taking sample 
points closer together than the theoretical coherence length 
for a given turbulence structure parameter, Cn?. Tilt 
problems in the FFT were suppressed by choosing the aperture 
array no larger than one-half the size of the FFT array, with 
one-fourth size giving the best results. 

Accuracy of the simulation was verified at various 
checkpoints within the coding. Aperture diffraction patterns 
and autocorrelations were compared to expected results. 
Turbulence effected diffraction patterns and atmospheric 
mutual coherence functions were presented to show trends for 
increasing turbulence. Calculated coherence length values 
were approximately 30% larger than theoretical coherence 
values when using a 512 by 512 FFT array with a 128 by 128 
aperture array. These results were obtained with only a 
Gaussian distributed random phase screen in the simulation. 

The addition of a Gaussian distributed random amplitude 


screen (for the complete Rytov approximation) to the 
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simulation brought calculated coherence values to within 3% 
of theoretical values. This showed that a numerical error was 
not present in the coding, but that turbulence effects were 
underestimated without the random amplitude term included. 
A single step, Fraunhofer algorithm must include the amplitude 
term to simulate turbulence effects on an optical wave 
correctly. Unity saturation of the normalized intensity 
variance occurred for Ca? values of approximately 10-11. 

Coherence length values used two types of filtering. The 
classical method and the Martin and Flatté method gave 
comparable results. The computational efficiency of the 
Martin and Flatté is significant. 

Future endeavors with this simulation should include the 
incorporation of Fresnel propagation into the coding. Also, 
further studies with the complete Rytov approximation in the 
simulation code should be pursued to obtain better statistical 


results. 
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APPENDIX A. RANDOM NUMBER GENERATOR SUBROUTINE 


=o om a ow om oe ow SF ee ow ae a=‏ سی سی سی سی سی am am ae oe‏ مت TU O a O O a‏ کہ ae a ae om ae ee oe CHAUD ow a eo ewe ao‏ مس سر سور سر سب سب می سی می روي a a‏ سي سب سر مس me eee ae a a aw a‏ سے سے 


subroutine RAN(iran,r) 


This is the random number generator algorithm from Dr. 
Harrison's notes, (Ref. 7). 


Variables: 
iran...input seed value (5 digit integer) 
CS ae returned random number uniformly distributed from 0 to 1 


iran=iran*99947 
r=0.5 + real(iran)*2.328306e-10 


return 
end 
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APPENDIX B. GAUSSIAN DISTRIBUTION SUBROUTINE 


e Te ae و سه‎ ee ee A A ee ee EA کچ‎ A Cr ER et ue A Er Cr Ce Cp Ce CE CD Cr ER GU A € E» UR UP GU auum E E CP D GP سي ټس سب‎ A MP CD A CEP CRF OEP GE aa ame a a o me o Mur می‎ ee جرد‎ 


subroutine GAUSS (narray, iran) 


This subroutine changes uniformly distributed random numbers, 
provided by subroutine RAN, to Gaussian distributed random 
numbers with unity variance. This subroutine is found in 
Knuth, [Ref. 8]. 
Variables: 

phaser..... real phase screen 2-D array 

phasel..... imaginary phase screen 2-D array 

narray..... dimension of the 2-D arrays 

7۶90-۹ input seed value for RAN (5 digit integer) 

E uniformly distributed random numbers 


xl & x2....Gaussian distributed random numbers 
common /blk2/phaser(512,512),phasei(512,512) 


do 10 i-1,narray 
do 10 j=1,narray 

20 call RAN(iran,r) 
vi=2.*r-1. 
call RAN(iran,r) 
v2=2.*r-1. 
s=v1*v1 + v2*v2 
if (s.ge.1.0) goto 20 
scale=sqrt (-2.*alog(s)/s) 
xl=vl*scale 
x2=v2*scale 
phaser(i,j)=x1 
phasei (1,3) =x2 

10 continue 


return 
end 
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APPENDIX C. SIMULATION CODE 


This code simulates the propagation of a monochromatic optical 
wave in a turbulent atmosphere.  Two-dimensional FFT routines 
are used extensively for calculations. Only Fraunhofer 
propagation is present in the coding. Filtering is accomplished 
by a choice of two methods: (1) classical and (2) the Martin 

and Flatté method. 


Subroutines: 

INIT.....Sets electric field arrays to zero. A value of 0 
in the call initializes real and imaginary arrays. 

A value of 1 initializes only the imaginary array. 

SQUARE...Establishes the planar electric field for the case 
of a square aperture. 

CIRCLE...Establishes the planar electric field for the case 
of a circular aperture. 

PLOT. ee Gives a screen plot of the array chosen in the call 
using EGA graphics. Different colors are chosen to 
represent different orders of magnitude values. Most 
calls within this subroutine are NDP Fortran-386 
specific. 

XFORM....Takes the 2-D FFT of the two arrays given in the call. 
Returns values in the same arrays. Makes calls to 
subroutine FFT (1-D FFT routine) several times. A 
value of -1 in the call is for a direct transforn. 

À value of *1 in the call is for an inverse transform. 

FFT......1-D FFT routine supplied by Dr. Don Walters. This 
routine is used by subroutine XFORM multiple times 
to accomplish the 2-D transform. 

This is Dr. Walters FFT 

MAG..... .Calculates the electric field magnitude and the 
intensity field from the real and imaginary electric 
field arrays. 

MCFPLOT..Calculates the Mutual Coherence Function (MCF) from 
the intensity field and then plots it on the screen 
using EGA graphics. Many calls within this 
subroutine are NDP Fortran-386 specific. 

GAUSS....Creates the phase screens from Gaussian distributed 
random numbers with unit variance. Makes calls 
to subroutine RAN for uniformly distributed random 
numbers then transforms them to a Gaussian distribution. 
Elements of this subroutine provided by Knuth. (Ref. 8] 

RAN. ..... Generates uniformly distributed random numbers. This 
algorithm from Harrison [Ref. 7) 
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FILTER...Filters the phase screens according to the Kolmogorov 
power law. A shuffling of array values is necessary 
for proper filtering. 

MESH.....Meshes the random phase screen with the electric 
field. 


Variables: 
fieldr...Real 2-D electric field array. 
fieldi...Imaginary 2-D electric field array. 
fieldro..Real 2-D planar electric field at aperture. 
phaser...Real 2-D phase screen array. 
phasei...Imaginary 2-D phase screen array. 
fieldm...2-D array that represents the electric field magnitude. 
fieldm2..2-D array that represents the intensity field. 
narray...Dimension of electric field array 
nfield...Dimension of aperture array 
lran.....Input seed for random number generator (5 digit integer). 


CHI. E Refractive index structure parameter. 
0 e ae Wavelength of monochromatic electromagnetic wave. 
delz..... Physical distance between aperture plane and image plane. 


width....Physical width of aperture plane. 

delx.....Physical distance between sample points in aperture 
array. 

rhonot...Coherence length calculated by interpolation between 
points of MCF curve. 

ren... Real 1-D array used by subroutine FFT. 

Tif ٦۷٦ Imaginary 1-D array used by subroutine FFT. 

fphaser..Real 2-D array used in subroutine FILTER that represents 
the phase screen in a more convenient form for 
filtering and viewing. 

fphasei..Imaginary 2-D array used in subroutine FILTER that 
represents the phase screen in a more convenient form 
for filtering and viewing. ۱ 


common /b1k1/fieldr(512,512) fieldro(512,512), fieldi (512,512) 
common /b1k2/phaser (512,512) , phasei (512,512) 
common /bl1k3/fieldm(512,512) ,fieldm2 (512,512) 


Initialize arrays 
1 call INIT(narray,0) 
Input section 


write(*,*) ‘Enter dimension of array that required. 
read(*,*) narray 
write(*,*) ‘Enter dimension of planar electric field.’ 
write(*,*) ‘(Pixel width of the aperture) ' 
read(*,*) nfield 

2 write(*,*) "Choose aperture shape: 1) SQUARE' 
write(*,*) ' 2) CIRCLE' 
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read(*,*) ichoice 

if (ichoice .1t. 1 .or. ichoice .gt. 2) then 
write(*,*) ‘Try again!' 
goto 2 

endif 

write(*,*) ‘Enter the seed for the random number generator. ' 

write(*,*) ‘(Value must be a five digit integer) ' 

read(*,*) iran 

write(*,*) ‘Enter the value of Cn squared.’ 

read(*,*) cn2 

write(*,*) ‘Enter the wavelength of light." 

read(*,*) wvl 

write(*,*) "Enter the distance from aperture to observer.’ 

read(*,*) delz 

0461221000 . 

Write(*,*) 'Enter the width of the aperture in meters.' 

read(*,*) width 

delx=width/real (nfield) 

Load arrays with desired planar electric field 


if (ichoice .eq. 1) call SQUARE (narray,nfield) 
if (ichoice .eq. 2) call CIRCLE (narray,nfield) 


Plot planar electric field 
call PLOT (fieldr,narray,1) 
Take Fourier transform of the field 
call XFORM(fieldr,fieldi,narray,delx,-1.) 
Calculate magnitude of transformed field then plot it 
call MAG(narray) 
call PLOT(fieldm,narray,1) 
Set imaginary portion of field to zero 
call INIT (narray,1) 
Take Fourier transform of field intensity 
call XFORM(fieldm2,fieldi,narray,delx,*1.) 
Calculate and plot MCF of aperture 
call MCFPLOT (narray,delx,nfield,0) 
Load phase screen arrays with Gaussian random numbers 
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call GAUSS (narray,iran) 
Transform phase screen to frequency space 
call XFORM (phaser, phasei,narray,sqrt (real (narray)),-1.) 
call PLOT(phaser,narray,1) 
call PLOT(phasei,narray,1) 
Filter phase screen using Kolmogorov spectrum idea 
call FILTER (narray,cn2,wvl,delx,delz) 
call PLOT(phaser,narray,1) 
call PLOT(phasei,narray,1) 
Transform phase screen back to real space 
call XFORM(phaser,phasei,narray,sqrt(real(narray)*delx),*1.) 
call PLOT(phaser,narray,1) 
call PLOT(phasei,narray,1) 


Mesh phase screen with planar electric field 


call INIT (narray,1) 
call MESH (narray) 


call PLOT(fieldr,narray,1) 

Take fourier transform of this electric field 
call XFORM(fieldr,fieldi,narray,delx,-1.) 
call MAG (narray) 
call PLOT(fieldm,narray,1) 

Reset imaginary portion of field to zero 
call INIT(narray,1) 

Take Fourier transform to get diffraction pattern w/ turbulence. 
call XFORM(fieldm2,fieldi,narray,delx,+1.) 

Calculate atmospheric MCF and plot it. 
call MCFPLOT (narray,delx,nfield,1) 
goto 1 


stop 
end 
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subroutine INIT(narray,k) 


common /blki/fieldr(512,512),fieldro(512,512),fie1di(512,512) 
do 10 i=1,narray 
do 10 j=1,narray 
if (k .eq. 0) then 
fieldr(1,3)=0.0 
fieldro(i,j)=0.0 
endif 
fieldi(i,j)=0.0 
10 continue 


return 
end 
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subroutine MAG(narray) 


common /b1k1/fieldr(512,512) fieldro(512,512), fieldi (512,512) 
common /b1k3/fieldm(512,512) fieldm2(512,512) 


do 10 i=1,narray 
do 10 j=1,narray 
fieldm2(i,j)sfieldr(i,j)**2 * fieldi(i,j)**2 
fieldnm(i,j)ssqrt(fieldn2(i,j)) 
10 continue 


return 
end 


subroutine SQUARE(narray,nfield) 
common /blki/fieldr (512,512) ,fieldro(512,512) , fieldi (512,512) 


do 10 i=narray/2+1-nfield/2,narray/2+nfield/2 
do 10 j=narray/2+1-nfield/2,narray/2+nfield/2 
fieldr(i,j)=1.0 
fieldro(i,3)=1.0 
10 continue 


return 
end 
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subroutine CIRCLE(narray,nfield) 


common /blk1/fieldr (512,512) ,fieldro(512,512) , fieldi (512,512) 
narray2=narray/2 
nfield2=nfield/2 


do 10 i=narray2+1-nfield2,narray2+nfield2 
do 20 j=narray2+1-nfield2,narray2+nfield2 
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radius2= (real (i)-real (narray2))**2 + 
* (real (j)-real (narray2))**2 
if (radius2 .1t. real(nfield2)**2) then 
fieldr(i,j)=1.0 
fieldro(i,j)=1.0 
endif 
20 continue 
10 continue 


return 
end 


subroutine PLOT(field,narray,k) 


dimension field(512,512) 
dimension ndex(20) 
data ndex/4,4,12,12,14,14,10,10,2,2,3,3,9,9,1,1,8,8,0,0/ 


fmax=0.0 
do 100 i=1,narray 
do 100 j=1,narray 
x=field(i,j) 
if (x.gt.fmax) then 
fmax=x 
endif 
100 continue 


call set video mode(16) 
call ega set mode 4 


do 110 i=1,narray/k 
do 110 j=1,narray/k 

ix-i 

iy-j 

xmax-alogiO(field(i,j)/fmax) 

index=8*abs (xmax) +1 

if (index.ge.21) then 
icolor=0 

else 
icolor=ndex (index) 

endif 

call ega_put_pixel (ix,iy,icolor) 

110 continue 


call pause 
call set_video_mode (3) 


return 
end 


subroutine MCFPLOT(narray,delx,nfield,iaprture) 
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common /blk3/fieldm(512,512),fieldm2(512,512) 
dimension apmcf (513),ymcf (513) 

logical flag 

efold=0.36787944 

rhonot=0. 

flag=.true. 


do 10 i=1,513 
if (iaprture .eq. 0) then 
apmef (1)=0.0 
endif 
ymcf(i)=0.0 
10 continue 


fmcf=0.0 
do 20 i=1,narray 
do 20 j=1,narray 
xmcf=fieldm2(i,)3) 
1f (xmcf.gt.fmcf) then 
fmcf=xmcf 
endif 
20 continue 


do 30 146 

if (iaprture .eq. 0) then 
apmcf (1)=fieldm2(1,1)/fmcf 
ymcf (1)=apmcf (1) 

else 
ymcf (1)2fieldm2(1,i)/(apnmcf(i)*fmcf) 

endif 

30 continue 


do 40 i=i,nfield 
if (iaprture.eq.1) then 
if (flag) then 
if (ymcf(i).lt.efold) then 
rhonot=real (1-1) *delx 
* + (efold-ymcf(i))*delx/(ymcf (i-1)-ymcf(1)) 
flag=.false. 
endif 
endif 
endif 
40 continue 


call set, video, mode (16) 
call ega set mode 4 
call locate(20,1) 
if (iaprture .eq. 0) then 
call write string('Aperture MCF vs. Length') 
else 
call write string('Atmospheric MCF vs. Length") 
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endif 

call locate(1,1) 

call write string('1.0') 

call locate(25,24) 

call write string('length') 

call locate(37,24) 

call write string(nfield) 

call ega draw, line(40,20,40,320,15) 
call ega draw line(40,320,600,320,15) 


ix-40 

1y=20 

do 50 i=1,nfield 
1x2=40+560*1/nfield 
iy2=20+int(300.*(1.-ymcf (i+1))) 
call ega_draw_line(ix,iy,ix2,iy2,15) 
ix=1x2 
1y=1y2 

50 continue 


call pause 
call set_video_mode (3) 


write(*,*) 'The coherence length is ',rhonot,' meters.' 


return 
end 


subroutine XFORM(fieldr,fieldi,narray,fftnorm,sign) 


common /blk4/re(512),rim(512) 
dimension fieldr (512,512) ,fieldi (512,512) 
data re/512*0./,rim/512*0./ 


m=int (alog(real(narray))/alog(2.)) 


do 30 i=l,narray 
do 40 j=1,narray 
re(j)2fieldr(i,j) 
rim(j)=fieldi(i,3) 
40 continue 
call FFT(m,fftnorm,sign) 
do 50 j=1,narray 
fieldr(i,j)=re(j) 
fieldi(i,j)=rim(j) 
50 continue 
30 continue 
do 60 j=1,narray 
do 70 i-1,narray 
re (i)=fieldr(i,3) 
rim(i)=fieldi (i,j) 
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70 continue 
call FFT(m,fftnorm,sign) 
do 80 i=1,narray 
fieldr(i,3)=re(i) 
fieldi(i,j)=rim(i) 
80 continue 
60 continue 


return 
end 


subroutine FFT(m,fftnorm,sign) 


common /blk4/re(512),rim(512) 
pi=3.141592653589792*sign 


n=2**m 
n1=n-1 
j=1 


do 200 i=1,n1 
if (i.1t.j) then 
t=re(j) 
re(j)=re(i) 
re(i)=t 
t=rim(3) 
rim(3)=rim(i) 
rim(1)=t 
endif 
k=n / 2 
do 201 while (k.1t.j) 
j=j-k 
k=k/2 
201 continue 
j=j+k 
200 continue 
le=1 
do 202 8 
lei=le 
le=letle 
ure=1. 
uim=0. 
ang=pi/lel 
wre=cos (ang) 
Wim=sin (ang) 
do 203 j=1,lel 
do 204 66 
ip=itlel 
tre=re (ip) *ure-rim(ip) *uim 
tim=re (ip) *uimtrim (ip) *ure 
re(ip)=re(i)-tre 
rim(ip)=rim(i)-tim 
re(i)=re(i)+tre 
rim(i)=rim(i)+tim 
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204 continue 
t-ure*wre-uim*wim 
uim-ure*wimtuim*wre 
ure-t 

203 continue 

202 continue 

if (sign.gt.0.0) then 
pts=1.0/real (n*fftnorm) 


else 
pts=fftnorm 

endif 

do 205 i=1,n 


re(i)=re(i)*pts 
rim(i)=rim(i)*pts 
205 continue 


return 
end 


ست ee‏ ہے ہجے ce ee ee‏ سس CURED‏ سم سم سس eee ee ee TÍ‏ ټس ټس نټ جس TÍ‏ نم TÍ Que Ga ee ee ee ee ee ee ee es es ee ee A ee ee ee‏ دوس سک ee CU‏ سر مس E‏ سه O A E AO A‏ جہن مش مس مس سه مس بس کسی SO ds ds AS‏ سے 


subroutine GAUSS (narray,iran) 
common /b1k2/phaser (512,512) ,phasei(512,512) 


do 10 i=1,narray 
do 10 j=1,narray 

20 call RAN(iran,r) 
vl=2.*r-1. 
call RAN(iran,r) 
v2=2.*r-1. 
s=v1*v1 + v2*v2 
if (s.ge.1.0) goto 20 
scale=sqrt (-2.*alog(s)/s) 
x1=v1*scale 
x2=v2*scale 
phaser (i,j)=x1 
phasei (1,))=x2 

10 continue 


return 
end 


e‏ چجھے کہ اس ee a a ee ee‏ سم سر e‏ مس ee ee daa CS dis CA‏ جج کہ جنس جک وک ee ee ee ee a‏ مس مس مس جس A‏ مس لصي A‏ مس سے سے سے سے سے سس سس سس سس e le ci ce ci e ci‏ سه مس سے جس مس سے سے س ج س سے 


subroutine RAN(iran,r) 
kran=99947 


iran=iran*kran 
r=0.5 + real(iran)*2.328306e-10 


return 
end 


سے سم مس ہتے کس ee eo em om om om om om om om om om om‏ ہت ome ome ome a am a aD am am am am am am a am am a an dam a‏ مسر مت مس سر مس ہت ہت ہت مس ae ae ® a a a em GENER‏ مس مس سے سے یت سے سے سے سے سے سے 


subroutine FILTER (narray,cn2,wvl,delx,delz) 
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common /b1kx2/phaser (512,512) phasei (512,512) 
dimension fphaser (512,512) fphasei (512,512) 
dimension aaa(512,512) 

pi=3.141592653589792 

tpi=2.*pi 

narray2=narray/2 

npivot=narray2+1 

power=-11./6. 

factor=sqrt (1.303*(tpi**3)*cn2*delz)/wvl 
factor2=(tpi/real(delx*narray))**power 


do 10 i=1,narray 
do 10 j=1,narray 
fphaser (1,3)=0. 
fphasei(i,j)=0. 
10 continue 


do 20 i=l,narray 
do 20 j=1,narray 
if (i.le.npivot) then 
if (j.le.npivot) then 
fphaser(i-1*narray2,j-1*narray2)-phaser (i, j) 
fphasei (i-l+narray2,)-l+narray2)=phasei(i,J) 
else 
fphaser (i-1+narray2, j-1-narray2)=phaser (i,j) 
f phasei (i-l+narray2,j-1-narray2)=phasei (i,j) 
endif 
else 
if (j.le.npivot) then 
f phaser (i-1-narray2,j-ltnarray2) =phaser (i,j) 
fphasei(i-1-narray2,j-1*narray2)-phasei(i,]j) 
else 
fphaser (i-1-narray2,j-1-narray2)-phaser(i,j) 
f phasei (i-1-narray2,j-1-narray2)=phasei (i,j) 
endif 
endif 
20 continue 


do 30 i=1,narray 
do 30 j=1,narray 

freq=sqrt (real (i-narray2)**2 + real(j-narray2) **2) 

if (i.eq.narray2.and.j.eq.narray2) then 
fphaser(i,j)=0. 
fphasei(i,))=0. 

else 
fphaser (i,))=fphaser(i,))*factor*factor2*(freq)**power 
fphasei(i,))=fphasei(i,))*factor*factor2*(freq)**power 

endif 

30 continue 
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do 40 i=1,narray 
do 40 j=1,narray 
if (i.1t.narray2) then 
if (j.1t.narray2) then 
phaser (i+1+narray2,j+1+narray2)=fphaser (i,j) 
phasei(iti+tnarray2,j+titnarray2)=fphasei(i,j) 
else 
phaser (iti*narray2,j*1-narray2)-fphaser (i, j) 
phaseil (i+1+narray2,j+1-narray2)=fphasei (i,)3) 
endif 
else 
if (j.1t.narray2) then 
phaser (i+1-narray2,j+1+narray2)=fphaser (1,3) 
phasei(i*i-narray2,j*1*narray2)-fphasei (i,j) 
else 
phaser (it1-narray2, j+1-narray2) =fphaser (i,j) 
phasei(i*i-narray2,j*1-narray2)-fphasei(i,j) 
endif 
endif 
40 continue 


return 
end 


subroutine MESH(narray) 


common /blk1/fieldr (512,512),fieldro(512,512),fieldi(512,512) 
common /blk2/phaser(512,512),phasei(512,512) 


do 10 i=1,narray 
do 10 j=1,narray 
ercosphisfieldro(i,j)*cos(phaser (i, j))*exp(phasei(i,j)) 
eicosphi=fieldi(i,j)*cos(phaser(i,j))*exp(phasei(i,j)) 
ersinphi=fieldro(i,j)*sin(phaser(i,j))*exp(phasei(i,j)) 
eisinphi=fieldi(i,j])*sin(phaser(i,j))*exp(phasei(i,j)) 
fieldr(i,j)sercosphi-eisinphi 
fieldi(i,j)sersinphi*eicosphi 
10 continue 


return 
end 
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